
#include "../Lattice.h"
#include <string>
#include <sstream>

using namespace std;

#define Nx				512
#define Ny				512

#define Re				10000						// Reynolds
#define RMS				0.04						//RMS de la velocidad inicial LBM (RMS de PS es 1). se normaliza con este valor.
#define Nu				(RMS*Nx/(2*M_PI*Re))		// viscosidad (unidades LBM)
#define Tau				((6.0*Nu+1)/2.0)			//0.500977848 en 512 //0.500488924087 en 256 //podriamos hacer que se calcule
#define Omega			(1./Tau)					//factor de viscosidad de LBM
#define Rho				1.0							//densidad cte

#define LBM2PStime		(2*M_PI*RMS/Nx)				//factor de transformacion de tiempo LBM a PS
#define Ntimes			10000						//(10000*(1./1024.)/LBM2PStime) // pasos temporales. nstep(PS)*dt(PS)/LBM2PStime
#define Nglobals		50*((1./1024.)/(LBM2PStime))//periodo de escritura de globals. nstg(PS)*dt(PS)/LBM2PStime
#define Ngraph			(Nglobals)					//periodo de escritura de graficos (i.e. vorticidad).

#define Delta			5.0							//dx del shear layer
#define noiseVolume 	0.002						//0.004 //ruido aleatorio (conviene que produzca un cambio del 10% de la energia)
#define noiseKMax		60							//maximo k del ruido

const string strV = "vorticity.";
const string strS = "spectrum.";

void print_globals(Lattice lat) {
	cout << "RMS real velocity = " << lat.RMSvelocity() << endl;
	cout << "energy = " << lat.energy()/(Nx*Ny) << endl;
	cout << "energy_k = " << lat.energy(FOURIER_SPACE) << endl;
	cout << "enstrophy = " << lat.enstrophy() << endl;
	cout << "palinstrophy = " << lat.palinstrophy() << endl;
}

void writeParameters(ostream& outStream) {
	outStream << "&INPUTS NX = " << Nx << ", NY = " << Ny << ", OMEGA = " << Omega << " /" << endl << endl;
	outStream << "Re = " << Re << endl;
	outStream << "Nu = " << Nu << endl;
	outStream << "Tau = " << Tau << endl;
	outStream << "LBM2PStime = " << LBM2PStime << endl;
	outStream << "Ntimes = " << Ntimes << endl;
	outStream << "Ngraph = " << Ngraph << endl;
	outStream << "Delta = " << Delta << endl;
	outStream << "Noise = " << noiseVolume << endl;
}

template <class T> inline std::string toString (const T& t) {
	std::stringstream ss;
	ss << t;
	return ss.str();
}

int main(int argc, char *argv[]) {
	Lattice lattice;
	lattice.init(Nx, Ny, true, true, Rho, Omega, 0);
	cout.precision(8);

	for(int i=0; i < Nx; i++)
		for(int j=0; j < Ny; j++) {
			// double a = (i < Nx/2.)*1.0 + (i >= Nx/2.)*(-1.0);
			// double b = (j < Ny/2.)*1.0 + (j >= Ny/2.)*(-1.0);
			//vec<double> _u(1.0, 0.0); //le agregué el 1.0 + para que la velocidad vaya entre 1 y -1
			vec<double> _u(1.0 + 1.0*(tanh((j-3.0*Ny/4.0)/(Delta))-tanh((j-Ny/4.0)/(Delta))), 0.0); //le agregué el 1.0 + para que la velocidad vaya entre 1 y -1
			lattice.setU(i, j, _u, 1.0);
		}

	//print parameters
	cout << "--- Parameters ---" << endl;
	writeParameters(cout);
	ofstream file("parameters.lbm");
	writeParameters(file);
	file.close();
	cout << endl << "parameters.lbm written." << endl << endl;
	cout << "--- Globals ---" << endl;
	print_globals(lattice);

	//noise
	double energyOld = lattice.energy();
	lattice.introRandomFourierNoise(noiseVolume, noiseKMax);
	double energyNew = lattice.energy();
	cout << endl << "--- Noise added ---" << endl;
	print_globals(lattice);
	cout << (energyNew-energyOld)*100./energyOld << "% of energy added" << endl;

	//PS normalization
	lattice.normalizeLattice(1.0/lattice.RMSvelocity());
	cout << endl << "--- Lattice normalized for PS ---" << endl;
	print_globals(lattice);

	//write PS input
	lattice.writeVorticity("kvorticity.lbm", FOURIER_SPACE);
	cout << endl << "kvorticity.lbm written." << endl;

	//LBM normalization
	lattice.normalizeLattice(RMS/lattice.RMSvelocity());
	cout << endl << "--- Lattice normalized for LBM ---" << endl;
	print_globals(lattice);

	//LBM integration
	//lattice.relaxRho();
	cout << endl << "--- Starting LBM ---" << endl;
	for(int n=0; n<Ntimes; n++) {
		if (n % int(Nglobals) == 0)
			cout << n*LBM2PStime << "\t\t" << lattice.getPopulation() << "\t\t" << lattice.energy()/(Nx*Ny) << "\t\t" << lattice.enstrophy() << "\t\t" << lattice.palinstrophy() << endl; //escribo el tiempo en unidades de PS
		/*if (n % int(Ngraph) == 0) {
			lattice.writeVorticity(strV.append(toString(n)), REAL_SPACE);
			lattice.writeSpectrum(strS.append(toString(n)));
		}*/
		lattice.update();
	}
	return 0;
}
